function v = f31(theta,phi,sigma,T,X)
v=0;
x=[cos(phi)*sin(theta); sin(phi)*sin(theta); cos(theta)];
for j=1:T
    d=X(:,j)'*x;
    v=v+exp(-(1-d^2)/(2*sigma*sigma));
end

